Especialização em Estatística Aplicada (PPGEAD): Módulo de Estatística Espacial
Fonte: Storymap UFRRJ
Introdução à Análise Estatística Espacial
O que é Análise Estatística Espacial ?
São métodos estatísticos que levam em consideração a localização espacial do fenômeno estudado;
Segundo Bailey & Gatrell (1995), “Análise estatística espacial pode ser realizada quando os dados são espacialmente localizados e se considera explicitamente a possível importância de seu arranjo espacial na análise ou interpretação dos resultados”
A primeira pergunta a ser feita é: A distribuição dos dados apresenta um padrão aleatório ou apresenta uma agregação definida (clusters) ?
Origem da Estatística Espacial
- Dr. John Snow (1813-1858) \(\rightarrow\) Considerado pai da Epidemiologia Moderna
Mapeamento dos casos de coléra (\(\bullet\)) e as bombas de água (X) em Londres, 1854.
Objetivos da Estatística Espacial
Análise de padrões espaciais e espaço-temporais (ex: Análise Exploratória Espacial, Correlação Espacial)
Modelagem Espacial (ex: Modelos Estatísticos de Regressão Espacial)
Dependência Espacial ou Autocorrelação Espacial
“Independência é um pressuposto muito conveniente que faz grande parte da teoria estatística matemática tratável. Entretanto, modelos que envolvem dependência estatística são frequentemente mais realísticos. . . . Neste tipos de dados a dependência está presente em todas as direções e fica mais fraca a medida em que aumenta a dispersão na localização dos dados.” (Cressie N.,1991)
“Todas as coisas se parecem, coisas mais próximas são mais parecidas que aquelas mais distantes.” (Tobler, 1979). Também conhecida como \(1^a\) Lei da Geografia
Tipologia dos Dados Espaciais
Os diferentes tipos de dados espaciais são tradicionalmente classificados de acordo com uma tipologia. Esta caracterização diz respeito a natureza estocástica da observação.
Noel Cressie (1993) divide a estatı́stica espacial em 3 grandes áreas:
Dados de processos pontuais
Dados de geoestatística
Dados de área
Existem métodos estatı́sticos diferentes para descrever ou analisar estes tipos de dados. Exemplo:
Padrões Pontuais
O principal interesse está no conjunto de coordenadas geográficas representando as localizações exatas de eventos.
Exemplos: Localização de crimes, localização da residência dos casos de dengue, localização de espécies vegetais, etc.
Neste caso, o dado aleatório de interesse é a localização espacial do evento.
Estimativas de Kernel (ou Mapas de Calor)
- Estimador de intensidade de distribuição de pontos:
\[\hat{\lambda}_{\tau}(u) = \dfrac{1}{\tau^2}\sum k(\dfrac{d(u_i , u)}{\tau}), d(u_i , u) \leq \tau\] *Fonte: Referência científica: Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds) “Análise Espacial de Dados Geográficos”. Brasília, EMBRAPA, 2004
- Localização da ocorrência de casos de dengue em Belo Horizonte/MG
Geoestatística
Podemos definir como sendo uma análise de um atributo espacialmente contínuo amostrado em localizações fixas.
Os dados compreendem um conjunto de localizações (em geral latitudes e longitudes), mas atribuidos a eles uma medida, como por exemplo: volume de chuva, concentração de poluentes no ar, número de ovos de Aedes postado em ovitrampas, etc.
A determinação experimental do semivariograma, para cada valor de \(x_i\), considera todos os pares de amostras \(z(x)\) e \(z(x+h)\), separadas pelo vetor distância \(h\), a partir da equação:
\[\hat{\gamma}(h) = \dfrac{1}{2N(h)}\sum^{N(h)}_{i=1}[z(x_i) - z(x_i + h)]^2 \]
- Sendo \(\hat{\gamma}(h)\) o semivariograma estimado e \(N(h)\) o número de pares de valores medidos, \(z(x)\) e \(z(x+h)\), separados pelo vetor \(h\).
- Esta fórmula, entretanto, não é robusta. Podem existir situações em que variabilidade local não é constante e se modifica ao longo da área de estudo (heteroscedasticidade).
- Exemplo: Mapa sobre o teor de argila no solo.
Dados de Área
Na análise de áreas o atributo estudado é em geral resultado de uma medida sumáraria (ex: contagem, taxas, médias, etc.) em uma determinada área bem definida, por exemplo um polígono.
Tal polígono cuja forma pode ser complexa bem como as relações de vizinhança.
O objetivo é a detecção e explicação de padrões e tendências observados entre áreas.
Mapa Temático
- Taxas de câncer de pulmão na população branca masculina nos Estados Unidos, por condados no ano de 1998
Matriz de Vizinhança
Moran Global, Moran Local e Lisa Map
Moran Global
Moran Local
- Desigualdade no nível distrital na cobertura de saúde reprodutiva, materna, neonatal e infantil na Índia
Modelos de Regressão Espacial
Hipótese de independência das observações em geral é Falsa \(\rightarrow\) Dependência Espacial
Efeitos Espaciais \(\rightarrow\) Se existir forte tendência ou correlação espacial, os resultados serão influenciados, apresentando associação estatística onde não existe (e vice-versa);
Como verificar ? \(\rightarrow\) Medir a autocorrelação espacial dos resíduos da regressão (Índice de Moran dos resíduos).
Autocorrelação espacial constatada ! E agora ? \(\rightarrow\) Utilizar Modelos de regressão que incorporam efeitos espaciais:
- Modelos Globais: utilizam um único parâmetro para capturar a estrutura de correlação espacial \(\rightarrow\) ex: Spatial Error Models (CAR)
\[y_i = \beta_0 + \sum^{p}_{k} \beta_k x_{ik} + \varepsilon_i \text{ , sendo } \varepsilon_i = \lambda W + \xi\]
Os efeitos da autocorrelação espacial são associados ao termo de erro. Sendo \(W\) a matriz de vizinhaça, \(\varepsilon\) a componente do erro com efeitos espaciais, \(λ\) é o coeficiente autoregressivo e \(\xi\) é a componente do erro com variância constante e não correlacionada.
A hipótese nula para a não existência de autocorrelação é que \(\lambda = 0\), ou seja, o termo de erro não é espacialmente correlacionado.
- Modelos Locais: parâmetros variam continuamente no espaço ex: Geographically Weighted Regression (GWR)
\[y_i = \beta_{0}(u_i, v_i) + \sum^{p}_{k} \beta_k (u_i, v_i) x_{ik} + \varepsilon_i \text{ , sendo } (u_i, v_i) \text{ as coordenadas geográficas}\] Fonte: ARDILLY, Pascal et al. Manuel d’analyse spatiale. 2018.
Como trabalhar com Análise Estatística Espacial: Algumas ferramentas
Serão feitas cinco aplicações da Estatística Espacial utilização o pacote estatístico R:
Clique aqui para baixar os dados das aplicações
Aplicação I: Utilizando a biblioteca tmap para construção de mapas temáticos
library(tmap)
data(World)
tmap_style("classic")
# Desenhando um rápido mapa temático mundial para esperança de vida.
qtm(World, fill = "life_exp")Aplicação II: Baixando e construindo mapas a partir da biblioteca geobr
Bibliotecas que serão utilizadas:
library(ggplot2)
library(dplyr)
library(viridis)
library(geobr)
library(sf)
library(maptools)
library(leaflet)
library(ggspatial)- Para acessar os dados dos limites territoriais de todos os estados brasileiros é necessário utilizar a função read_state.
Primeiro, vamos fazer um gráfico apenas com as geometrias.
- Para a construção de um mapa onde cada estado recebe uma cor de acordo com a sua região geogrática, procedemos da seguinte forma:
- Agora utilizaremos dados relativos a porcentagem de municípios com rede de esgoto de acordo com a unidade da federação (fonte dos dados ) para fazer um gráfico. Vamos associar esses dados a tabela de acordo com a variável State e padronizaremos a porcentagem para variar de 0 a 100.
# Entrando com os dados observados na wikipedia
acesso_san <- data.frame(code_state = c(12, 27, 16, 13, 29, 23, 53, 32, 52, 21, 51,
50, 31, 15, 25, 41, 26, 22, 33, 24, 43, 11, 14, 42, 35, 28, 17), com_rede = c(0.273,
0.412, 0.313, 0.177, 0.513, 0.696, 1, 0.974, 0.28, 0.065, 0.191, 0.449, 0.916,
0.063, 0.731, 0.421, 0.881, 0.045, 0.924, 0.353, 0.405, 0.096, 0.4, 0.352, 0.998,
0.347, 0.129))# Construindo o mapa com ggplot
brasil.ufs |>
left_join(acesso_san, by = "code_state") |>
mutate(com_rede = 100 * com_rede) |>
ggplot(aes(fill = com_rede), color = "black") + geom_sf() + scale_fill_viridis(name = "Municípios com rede de esgoto (%)",
direction = -1) + xlab("Longitude") + ylab("Latitude") + annotation_scale(location = "bl") +
annotation_north_arrow(location = "br") + theme_minimal()- Uma forma alteranativa de apresentar esses mesmos dados se dá pela apresentação de círculos com raios proporcionais a porcentgem de municípios com rede de esgoto no centroide de cada geometria (nesse caso, UF).
# Juntando o banco com os dados + malha
coord_pontos <- brasil.ufs |>
left_join(acesso_san, by = "code_state") |>
mutate(com_rede = 100 * com_rede) |>
st_centroid()# Construindo o mapa com ggplot
ggplot(brasil.ufs) + geom_sf() + geom_sf(data = coord_pontos, aes(size = com_rede),
col = "blue", alpha = 0.65, show.legend = "point") + scale_size_continuous(name = "Municípios com rede de esgoto (%)") +
xlab("Longitude") + ylab("Latitude") + annotation_scale(location = "bl") + annotation_north_arrow(location = "br") +
theme_minimal()- Uma alternativa interativa para trabalhar com mapas é com a utilização do pacote leaflet
data.frame(st_coordinates(coord_pontos), com_rede = coord_pontos$com_rede, UF = coord_pontos$name_state) |>
leaflet() |>
addTiles() |>
addCircleMarkers(~X, ~Y, label = ~as.character(paste0(UF, ": ", com_rede, "%")),
labelOptions = labelOptions(textsize = "13px"), radius = ~com_rede/10, fillOpacity = 0.5)Aplicação III: Dengue em Dourados/MS - Parte 1: Análise exploratória
Nesta apresentação serão utilizados os dados que fizeram parte da dissertação de Isis Rodrigues Reitman intitulada “SAÚDE E AMBIENTE URBANO: A RELAÇÃO DE INCIDÊNCIA DE DENGUE E AS DISPARIDADES ESPACIAIS EM DOURADOS – MS”, apresentada no Programa de Pós-Graduação em Geografia da Universidade Federal da Grande Douradosos/MS, em abril de 2016.
Essa aplicação também se encontra no Curso de Estudos Ecológicos ministrado para o curso de Pós-Graduação em Epidemiologia em Saúde Pública em 2019, pelos pesquisadores Oswaldo Gonçalves Cruz (PROCC/FIOCRUZ) e Wagner Tassinari (DEMAT/ICE/UFRRJ)
OBS: Os dados e as malhas geográficas utilizadas nessa apresentação, estão disponíveis no seguinte endereço: (link)
Biliotecas do R que serão utilizadas
Lendo a tabela da população por setor censitário e os shapes files do contorno e por setor censitário de Dourados/MS
pop2010 <- read_csv("dados/dengue_dourados/pop2010.csv")
setor.sf <- read_sf("dados/dengue_dourados/Setor_UTM_SIRGAS.shp", crs = 31981)
contorno.sf <- read_sf("dados/dengue_dourados/contorno.shp", crs = 31981)OBS: Um aspecto muito importante dos dados espaciais é o sistema de referência de coordenadas (CRS) que é usado. Por exemplo, uma localização de (140, 12) não é significativa se você sabe onde está a origem e se a coordenada x está a 140 metros, quilômetros ou talvez graus de distância dela (na direção x). (link)
Fazendo um join com as tabelas com os setores censitários + população
Lendo e plotando os casos de dengue georreferenciados em Dourados/MS
casos <- read_csv("dados/dengue_dourados/dengue_dourados.csv")
casos.sf <- st_as_sf(casos, coords = c("X", "Y"), crs = 31981)ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = casos.sf,
color = "red", size = 1) + theme_void()Lendo e plotando os os pontos de coleta de lixo georreferenciados em Dourados/MS
lixo <- read_csv2("dados/dengue_dourados/lixo_dourados.csv")
lixo.sf <- st_as_sf(lixo, coords = c("X", "Y"), crs = 31981)ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = lixo.sf,
color = "blue", size = 1) + theme_void()- Como podemos observar, existem alguns pontos de coleta de lixo fora do contorno de Dourados/MS
Uma forma de ficarmos só com os pontos dentro do polígono é utilizando o comando st_intersection.
ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = lixo2.sf,
color = "blue", size = 1) + theme_void()Utilizando as informações dos casos (pontos) + do lixo (ponto) + população de cada setor censitário (mapa temático)
ggplot(setor.sf) + geom_sf(aes(fill = pop)) + geom_sf(data = casos.sf, color = "red",
size = 0.7, aes(colour = "Caso"), show.legend = "point") + geom_sf(data = lixo2.sf,
color = "salmon", size = 1, aes(colour = "Lixo"), show.legend = "point") + scale_fill_distiller(palette = "PuBu",
direction = 1) + scale_colour_manual(values = c(Caso = "red", Lixo = "slamon")) +
theme_minimal()
#### Construindo os buffers
- Agora serão construidos buffers de 500m de distância ao redor de cada ponto de coleta de lixo. Isso é interessante para verificar se os casos de dengue ocorrem em um raio de até 500m de distância dos pontos de coleta de lixo. Ou seja, a pergunta é, será que a distância dos pontos de coeta de lixo influenciam a ocorrência do caso de dengue ?
Buffers: São polígonos que contornam um objeto a uma determinada distância. Sua principal função é materializar os conceitos de “perto” e “longe”.
lixo_buffer <- st_buffer(lixo2.sf, 500)
ggplot(setor.sf) + geom_sf(aes(fill = pop)) + geom_sf(data = lixo_buffer, color = "gray",
fill = "transparent", size = 0.4) + geom_sf(data = casos.sf, color = "red", size = 0.7) +
scale_fill_distiller(palette = "PuBu", direction = 1) + scale_colour_manual(values = c(Caso = "red",
Lixo = "slamon")) + theme_minimal()Represntando os casos e o lixo de forma interativa
Convertendo o dado de pontos (padrão pontual) para dados de área
## conta casos por setor
casos.sf$contador <- 1
casos <- setor.sf |>
st_join(casos.sf) |>
filter(CLASSI_FIN == 1) %>% ## seleciona somente os casos confirmados
group_by(ID) |>
summarise(casos = sum(contador))
st_geometry(casos) <- NULL ## remove atributos de geometria
## numero de depositos de Lixo por setor
lixo.sf$contador <- 1
lixo <- setor.sf |>
st_join(lixo.sf) |>
group_by(ID) |>
summarise(lixo = sum(contador))
st_geometry(lixo) <- NULL ## remove atributos de geometria
# Inserindo as contagens dos casos e de pontos de coleta de lixo no atributo com a geometria.
setor.sf <- setor.sf |>
left_join(lixo,by = 'ID') |>
left_join(casos,by = 'ID') Calculando a taxa de incidência e plotando o mapa temático dos pontos de coleta de lixo por setor censitário
setor.sf$tx <- (setor.sf$casos/setor.sf$pop) * 1000
setor.sf$tx[is.na(setor.sf$tx)] <- 0 # Transformando os missings em zero
summary(setor.sf$tx) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 0.000 1.736 4.065 4.311 56.399
Plotando a distribuição da incidência em Dourados/MS
library(wesanderson)
pal <- wes_palette("Moonrise3", 20, type = "continuous")
ggplot(setor.sf) + geom_sf(aes(fill = tx), color = "black") + scale_fill_gradientn(colours = pal) +
ggtitle("Taxa de incidência de Dengue") + theme_void()Kernel por atributos
Vamos plotar o kernel por atributos referente a taxa de incidência de dengue em Dourados/MS.
Primeiramente é necessário dissolver os polígonos em formato sf para obter o contorno. Nesse caso queremos preservar o atributo AREA
Extraindo os centróides dos polígonos em Dourados/MS
centroides <- st_centroid(st_geometry(setor.sf))
# Transformando em os centróides em formato sp
centroides.sp <- as.data.frame(as_Spatial(centroides))
names(centroides.sp) <- c("X", "Y")
plot(centroides.sp)Colocando os pontos no formato sp
centroides.ppp <- ppp(centroides.sp$X, centroides.sp$Y, dourados.w)
plot(centroides.ppp, pch = 19, cex = 0.5)Fazendo o kernel por atributo da taxa de detecção
kernel.tx <- density(centroides.ppp, 500, weights = setor.sf$tx, scalekernel = TRUE)
plot(kernel.tx)Construindo a matriz de vizinhança para verificar a autocorrelação espacial
Neighbour list object:
Number of regions: 284
Number of nonzero links: 1726
Percentage nonzero weights: 2.139952
Average number of links: 6.077465
- Iremos precisar da coordenadas dos centróides
setor.sp <- as(setor.sf, "Spatial") # convertendo em formato sp
coord <- coordinates(setor.sp) # coordenadas dos centroidas dos poligonos de dourados
class(setor.sp)[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
Verificando a malha de conectividade da vizinhança de Dourados/MS
viz.sf <- as(nb2lines(viz, coords = coord), "sf")
viz.sf <- st_set_crs(viz.sf, st_crs(setor.sf))
# Plota o grafo de conectividade por contiguidade
mapa.viz <- ggplot(setor.sf) + geom_sf(fill = "salmon", color = "white") + geom_sf(data = viz.sf) +
theme_minimal() + ggtitle("Vizinhança por \n conectividade") + ylab("Latitude") +
xlab("Longitude")
mapa.vizObtendo a correlação da taxa de incidência de dengue Dourados/MS
Moran I test under randomisation
data: setor.sf$tx
weights: pesos.viz
Moran I statistic standard deviate = 15.724, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.524720303 -0.003533569 0.001128660
Plotando o correlograma
Spatial correlogram for setor.sf$tx
method: Moran's I
estimate expectation variance standard deviate Pr(I) two sided
1 (284) 0.52472030 -0.00353357 0.00112866 15.7239 < 2.2e-16
2 (284) 0.23776816 -0.00353357 0.00048540 10.9525 < 2.2e-16
3 (284) 0.09536593 -0.00353357 0.00028358 5.8729 4.282e-09
4 (284) -0.02063378 -0.00353357 0.00019736 -1.2172 0.22351
5 (284) -0.07589158 -0.00353357 0.00016732 -5.5939 2.221e-08
6 (284) -0.06448448 -0.00353357 0.00016165 -4.7940 1.635e-06
7 (284) -0.02708340 -0.00353357 0.00016725 -1.8210 0.06861
8 (284) -0.02744543 -0.00353357 0.00018328 -1.7662 0.07735
1 (284) ***
2 (284) ***
3 (284) ***
4 (284)
5 (284) ***
6 (284) ***
7 (284) .
8 (284) .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Mapeando os polígonos que tiveram os p-valores mais significativos no Moran Local.
Moran Local (Lisa Map) da taxa de incidência Dourados/MS
Local Morans I Stand. dev. (N) Pr. (N) Saddlepoint Pr. (Sad) Expectation
1 1 -0.009885292 -0.01575255 1.0125682 -0.03600714 0.9712767 -0.003533569
2 2 0.226981101 0.46511580 0.6418485 0.70056077 0.4835772 -0.003533569
3 3 -0.010910944 -0.01829621 1.0145974 -0.04041673 0.9677609 -0.003533569
4 4 0.484675519 1.31014141 0.1901480 1.43930439 0.1500643 -0.003533569
5 5 0.018648578 0.05952726 0.9525322 0.09384037 0.9252360 -0.003533569
Variance Skewness Kurtosis Minimum Maximum omega sad.r
1 1 0.1625854 -0.05147857 8.753481 -57.75455 56.75455 -0.0001369880 -0.01569409
2 2 0.2456263 -0.04188294 8.752890 -70.87400 69.87400 0.0028119325 0.44472643
3 3 0.1625854 -0.05147857 8.753481 -57.75455 56.75455 -0.0001590844 -0.01822753
4 4 0.1388594 -0.05570264 8.753780 -53.41199 52.41199 0.0066304485 1.08297547
5 5 0.1388594 -0.05570264 8.753780 -53.41199 52.41199 0.0005595065 0.05929966
sad.u
1 1 -0.01569909
2 2 0.49831660
3 3 -0.01823490
4 4 1.59298211
5 5 0.05942125
[ reached 'max' / getOption("max.print") -- omitted 5 rows ]
setor.sf$MoranLocal <- summary(resI)[, 1]
library(scales)
ggplot(setor.sf) + geom_sf(aes(fill = MoranLocal), color = "black") + scale_fill_gradientn(colours = c("blue",
"white", "red"), values = rescale(c(min(setor.sf$MoranLocal), 0, max(setor.sf$MoranLocal))),
guide = "colorbar") + ggtitle("Moran local") + theme_void()Aplicação IV: Dengue em Dourados/MS - Parte 2: Modelagem (Modelos Linear, CAR e GWR)
Ajustando o modelo de regressão linear simples.
# Transformando os missings em zero
setor.sf$lixo[is.na(setor.sf$lixo)] <- 0
dourados.lm <- lm(tx ~ lixo, data = setor.sf)
summary(dourados.lm)
Call:
lm(formula = tx ~ lixo, data = setor.sf)
Residuals:
Min 1Q Median 3Q Max
-4.510 -4.115 -2.286 0.237 51.889
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.5103 0.4889 9.225 <2e-16 ***
lixo -0.3955 0.1945 -2.033 0.043 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.365 on 282 degrees of freedom
Multiple R-squared: 0.01445, Adjusted R-squared: 0.01095
F-statistic: 4.134 on 1 and 282 DF, p-value: 0.04298
Checando os residuos para verificar a presença de autocorrelação.
Moran I test under randomisation
data: dourados.lm$lmresid
weights: pesos.viz
Moran I statistic standard deviate = 15.578, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.520064269 -0.003533569 0.001129669
Ajustando o modelo CAR (Spatial Error Model)
library(spatialreg)
dourados.car <- errorsarlm(tx ~ lixo, data = setor.sf, listw = pesos.viz)
summary(dourados.car)
Call:errorsarlm(formula = tx ~ lixo, data = setor.sf, listw = pesos.viz)
Residuals:
Min 1Q Median 3Q Max
-12.53213 -2.34032 -1.09686 0.82021 31.62992
Type: error
Coefficients: (asymptotic standard errors)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.07199 1.54998 2.6271 0.008611
lixo -0.25129 0.14814 -1.6963 0.089826
Lambda: 0.8063, LR test value: 167.98, p-value: < 2.22e-16
Asymptotic standard error: 0.041575
z-value: 19.394, p-value: < 2.22e-16
Wald statistic: 376.13, p-value: < 2.22e-16
Log likelihood: -885.0643 for error model
ML residual variance (sigma squared): 25.435, (sigma: 5.0433)
Number of observations: 284
Number of parameters estimated: 4
AIC: 1778.1, (AIC for lm: 1944.1)
Checando os residuos para verificar a presença de autocorrelação
Moran I test under randomisation
data: dourados.car$carresid
weights: pesos.viz
Moran I statistic standard deviate = 0.78357, p-value = 0.2166
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.023005067 -0.003533569 0.001147099
Ajustando o modelo GWR (Geographically Weighted Regression)
# Biblioteca para ajustar o modelos GWR
library(spgwr)
# Estimando a largura de banda “ideal” para o kernel
GWRbanda <- gwr.sel(tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X, centroides.sp$Y),
adapt = T)Adaptive q: 0.381966 CV score: 14699.46
Adaptive q: 0.618034 CV score: 15071.4
Adaptive q: 0.236068 CV score: 14204.85
Adaptive q: 0.145898 CV score: 13336.85
Adaptive q: 0.09016994 CV score: 12151.44
Adaptive q: 0.05572809 CV score: 10995.75
Adaptive q: 0.03444185 CV score: 10093.37
Adaptive q: 0.02128624 CV score: 9758.837
Adaptive q: 0.01315562 CV score: 9295.706
Adaptive q: 0.008130619 CV score: 8996.895
Adaptive q: 0.005024999 CV score: 8988.675
Adaptive q: 0.00638842 CV score: 9000.179
Adaptive q: 0.00310562 CV score: 9842.835
Adaptive q: 0.004291861 CV score: 9067.819
Adaptive q: 0.005545779 CV score: 8976.568
Adaptive q: 0.005867639 CV score: 8980.638
Adaptive q: 0.005586469 CV score: 8976.67
Adaptive q: 0.005505089 CV score: 8976.598
Adaptive q: 0.005545779 CV score: 8976.568
# Ajustando o modelo GWR
dourados.gwr = gwr(tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X, centroides.sp$Y),
adapt = GWRbanda, hatmatrix = TRUE, se.fit = TRUE)
dourados.gwrCall:
gwr(formula = tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X,
centroides.sp$Y), adapt = GWRbanda, hatmatrix = TRUE, se.fit = TRUE)
Kernel function: gwr.Gauss
Adaptive quantile: 0.005545779 (about 1 of 284 data points)
Summary of GWR coefficient estimates at data points:
Min. 1st Qu. Median 3rd Qu. Max. Global
X.Intercept. -0.324952 1.412188 3.099163 5.397652 37.581278 4.5103
lixo -18.321709 -1.222217 -0.400929 -0.061492 1.555744 -0.3955
Number of data points: 284
Effective number of parameters (residual: 2traceS - traceS'S): 134.4484
Effective degrees of freedom (residual: 2traceS - traceS'S): 149.5516
Sigma (residual: 2traceS - traceS'S): 5.447405
Effective number of parameters (model: traceS): 100.4864
Effective degrees of freedom (model: traceS): 183.5136
Sigma (model: traceS): 4.917575
Sigma (ML): 3.952993
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 1904.233
AIC (GWR p. 96, eq. 4.22): 1687.144
Residual sum of squares: 4437.827
Quasi-global R2: 0.7140881
# Colocando a saída do modelo dentro de um objeto dataframe.
results <- as.data.frame(dourados.gwr$SDF)
head(results) sum.w X.Intercept. lixo X.Intercept._se lixo_se gwr.e
1 4.630090 6.489949 -1.8954895 2.411594 1.718037 -0.02159391
2 3.934954 7.883872 -1.2811203 2.504116 1.630176 1.98454873
3 3.062602 7.267798 -1.5625806 2.780695 1.625888 -1.85906371
4 4.139484 9.791258 -1.1687314 2.376803 1.277302 6.80957181
5 5.541437 4.586352 -0.6839813 2.090863 1.185501 -2.63010618
pred pred.se localR2 X.Intercept._se_EDF lixo_se_EDF pred.se.1
1 2.698970 2.477284 0.6403656 2.671425 1.903141 2.744191
2 7.883872 2.504116 0.5043081 2.773914 1.805815 2.773914
3 5.705218 2.164469 0.5288378 3.080292 1.801064 2.397673
4 8.622527 1.794393 0.4101644 2.632885 1.414921 1.987725
5 3.902371 1.524868 0.6665163 2.316137 1.313229 1.689160
coord.x coord.y
1 731406.5 7541547
2 730845.4 7541481
3 730957.7 7541239
4 730487.2 7541437
5 729874.1 7541446
[ reached 'max' / getOption("max.print") -- omitted 1 rows ]
Verificando a distribuição dos coeficientes de regressão para a variável lixo
Verificando a distribuição dos localR2
Incorporando alguns parâmetros de saída do modelo na tabela setor.sf
Mapa dos coeficientes de regressão para a variável lixo
map.lixo <- ggplot(setor.sf) + geom_sf(aes(fill = coef.lixo), color = "black") +
scale_fill_gradientn(colours = pal) + ggtitle("Distribuição dos coef. var. lixo") +
theme_void()
map.lixoChecando os residuos para verificar a presença de autocorrelação para o modelo GWR.
# Calculando os resíduos para o modelo GWR
results$residuos <- setor.sf$tx - results$pred
moran.test(results$residuos, pesos.viz)
Moran I test under randomisation
data: results$residuos
weights: pesos.viz
Moran I statistic standard deviate = 1.4806, p-value = 0.06935
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.046816835 -0.003533569 0.001156432
Mapeando os coeficientes de regressão para a variável lixo por significancia através do teste de wald
# Calculando a estatística de wald
setor.sf$wald.teste <- abs(results$lixo/results$lixo_se)
# Dicotomizando a estatística de wald
setor.sf$wald.teste <- ifelse(setor.sf$wald.teste < 2, 0, 1)
map.wald <- ggplot(setor.sf) + geom_sf(aes(fill = factor(wald.teste)), color = "black") +
scale_fill_manual(values = c("white", "purple"), labels = c("< 2", ">= 2"), name = "Wald") +
ggtitle("Coef. lixo significativos") + theme_void()
library(gridExtra)
grid.arrange(map.lixo, map.wald, ncol = 2)Mapa dos coeficientes de determinação regionalizados (\(R^2\) local).
ggplot(setor.sf) + geom_sf(aes(fill = localR2), color = "black") + scale_fill_gradientn(colours = pal) +
ggtitle("R² local") + theme_void()Verificando a distribuição dos preditos.
histdens <- function(x, titulo = "") {
densi <- density(x)
xli <- range(densi$x)
yli <- range(densi$y)
hist(x, col = "red", probability = T, xlim = xli, ylim = yli, main = titulo)
lines(densi, lwd = 2)
abline(v = median(x), lwd = 4, col = 4, lty = 2)
}
par(mfrow = c(2, 2))
hist.tx <- histdens(setor.sf$tx, "Tx Bruta")
hist.lm <- histdens(dourados.lm$fitted.values, "Pred LM")
hist.car <- histdens(dourados.car$fitted.values, "Pred CAR")
hist.gwr <- histdens(results$pred, "Pred GWR")Mapeando os valores observados e preditos dos modelos ajustados
library(colorspace) #
setor.sf$brks <- cut(setor.sf$tx, include.lowest = TRUE, right = TRUE, breaks = c(-4,
0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))
tx.bruta <- ggplot(setor.sf) + geom_sf(aes(fill = brks), color = "black") + ggtitle("Taxa Bruta") +
scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30,
l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") +
theme_void()
setor.sf$brks.lm <- cut(dourados.lm$fitted.values, lowest = TRUE, right = TRUE, breaks = c(-4,
0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))
pred.lm <- ggplot(setor.sf) + geom_sf(aes(fill = brks.lm), color = "black") + ggtitle("Taxa Predita - LM") +
scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30,
l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") +
theme_void()
setor.sf$brks.car <- cut(dourados.car$fitted.values, lowest = TRUE, right = TRUE,
breaks = c(-4, 0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10",
"> 10"))
pred.car <- ggplot(setor.sf) + geom_sf(aes(fill = brks.car), color = "black") + ggtitle("Taxa Predita - CAR") +
scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30,
l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") +
theme_void()
setor.sf$brks.gwr <- cut(results$pred, lowest = TRUE, right = TRUE, breaks = c(-4,
0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))
pred.gwr <- ggplot(setor.sf) + geom_sf(aes(fill = brks.car), color = "black") + ggtitle("Taxa Predita - GWR") +
scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30,
l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") +
theme_void()
library(gridExtra)
grid.arrange(tx.bruta, pred.lm, pred.car, pred.gwr, ncol = 2)Comparando a distribuição dos resíduos dos modelos ajustados
library(vioplot)
vioplot(dourados.lm$residuals, dourados.car$residuals, results$residuos, names = c("LM",
"CAR", "GWR"), col = "orange")
title("Gráficos violinos da distribuição dos resíduos")Mapeando os resíduos dos modelos ajustados
library(colorspace) #
setor.sf$brks.res.lm <- cut(dourados.lm$residuals, include.lowest = TRUE, right = TRUE,
breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5",
"> 5"))
res.lm <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.lm), color = "black") +
ggtitle("Resíduos - LM") + scale_fill_discrete_sequential(palette = "Purple-Yellow",
c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75",
drop = FALSE, name = "Taxa") + theme_void()
setor.sf$brks.res.car <- cut(dourados.car$residuals, include.lowest = TRUE, right = TRUE,
breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5",
"> 5"))
res.car <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.car), color = "black") +
ggtitle("Resíduos - CAR") + scale_fill_discrete_sequential(palette = "Purple-Yellow",
c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75",
drop = FALSE, name = "Taxa") + theme_void()
setor.sf$brks.res.gwr <- cut(results$residuos, include.lowest = TRUE, right = TRUE,
breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5",
"> 5"))
res.gwr <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.gwr), color = "black") +
ggtitle("Resíduos - GWR") + scale_fill_discrete_sequential(palette = "Purple-Yellow",
c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75",
drop = FALSE, name = "Taxa") + theme_void()
library(gridExtra)
grid.arrange(res.lm, res.car, res.gwr, ncol = 2)Aplicação V: Geoestatística com dados de pluviosidade na cidade do Rio de Janeiro/RJ
Biliotecas do R que serão utilizadas
Importando a tabela com a chuva acumulada média de 7 dias das últimas 24hs e 96hs das estações pluviométricas da cidde do Rio de Janeiro
Análise Gráfica Descritiva
0% 20% 40% 60% 80% 100%
0.00 1.92 4.24 7.36 11.84 36.00
Transformar os dados em um objeto espacial do R
- x - Longitude
- y - Latitude
Análise Gráfica Descritiva com os dados espaciais
# Point plot
spplot(pluvio["acumulado_24h"], scales = list(draw = T), key.space = "right", colorkey = T)Modelando o variograma experimental (ou empírico)
- width - Distância média entre amostras ou distância dos lags
- cutoff - Máxima distância
variogram.emp = variogram(acumulado_24h ~ x + y, pluvio, width = 1000, cutoff = 20000)
variogram.emp np dist gamma dir.hor dir.ver id
1 2 1385.875 20.381155 0 0 var1
2 6 2543.894 8.139327 0 0 var1
3 8 3648.158 5.261490 0 0 var1
4 15 4583.562 13.707344 0 0 var1
5 22 5491.324 26.228076 0 0 var1
6 15 6531.432 52.574023 0 0 var1
7 12 7512.398 42.179960 0 0 var1
8 20 8516.447 65.283855 0 0 var1
9 15 9432.920 36.022491 0 0 var1
10 19 10450.957 36.569808 0 0 var1
11 13 11492.778 76.150012 0 0 var1
12 19 12568.517 39.354425 0 0 var1
[ reached 'max' / getOption("max.print") -- omitted 7 rows ]
Ajustando o semivariograma teórico
- Sill - Semivariância estrutural ou contribuição (ponto máximo que chega ao plato no eixo de y)
- Range - Alcance (ponto máximo em x)
- Nugget - Efeito pepita, quando (y, 0)
## 1) Modelo Linear
lin.fit = fit.variogram(variogram.emp, model = vgm(psill = 240, model = "Lin", range = 5000,
nugget = 10))
## 2) Modelo exponencial
exp.fit = fit.variogram(variogram.emp, model = vgm(psill = 240, model = "Exp", range = 5000,
nugget = 10))
## 3) Modelo gaussiano
gau.fit = fit.variogram(variogram.emp, model = vgm(psill = 240, model = "Gau", range = 5000,
nugget = 10))
## 3) Modelo wave
wav.fit = fit.variogram(variogram.emp, model = vgm(psill = 240, model = "Wav", range = 5000,
nugget = 10))Validação cruzada
- RMSE (root mean squared error): é a medida que calcula “a raiz quadrática média” dos erros entre valores observados (reais) e predições (hipóteses).
## 1) Modelo Linear
cv.lin <- krige.cv(acumulado_24h ~ x + y, locations = pluvio, model = lin.fit)
summary(cv.lin)Object of class SpatialPointsDataFrame
Coordinates:
min max
x 634918.7 687966
y 7450222.8 7475508
Is projected: NA
proj4string : [NA]
Number of points: 30
Data attributes:
var1.pred var1.var observed residual
Min. : 0.6757 Min. : 9.001 Min. : 0.00 Min. :-10.41265
1st Qu.: 4.3063 1st Qu.:11.956 1st Qu.: 2.20 1st Qu.: -2.84581
Median : 6.7569 Median :19.120 Median : 5.00 Median : -0.44588
Mean : 7.2750 Mean :28.832 Mean : 7.34 Mean : 0.06504
3rd Qu.:10.4716 3rd Qu.:43.382 3rd Qu.:10.30 3rd Qu.: 1.60176
Max. :15.6127 Max. :80.832 Max. :36.00 Max. : 29.36110
zscore fold
Min. :-2.0312467 Min. : 1.00
1st Qu.:-0.5694160 1st Qu.: 8.25
Median :-0.1237453 Median :15.50
Mean : 0.0004815 Mean :15.50
3rd Qu.: 0.3209767 3rd Qu.:22.75
Max. : 4.1304203 Max. :30.00
plot(cv.lin$var1.pred ~ cv.lin$observed, cex = 1.2, lwd = 2)
abline(0, 1, col = "lightgrey", lwd = 2)
lm_lin <- lm(cv.lin$var1.pred ~ cv.lin$observed)
abline(lm_lin, col = "red", lwd = 2)r2_lin = summary(lm_lin)$r.squared
rmse_lin = hydroGOF::rmse(cv.lin$var1.pred, cv.lin$observed)
## 2) Modelo exponencial
cv.exp <- krige.cv(acumulado_24h ~ x + y, locations = pluvio, model = exp.fit)
summary(cv.exp)Object of class SpatialPointsDataFrame
Coordinates:
min max
x 634918.7 687966
y 7450222.8 7475508
Is projected: NA
proj4string : [NA]
Number of points: 30
Data attributes:
var1.pred var1.var observed residual
Min. : 0.06051 Min. : 12.38 Min. : 0.00 Min. :-15.0221
1st Qu.: 3.53749 1st Qu.: 18.23 1st Qu.: 2.20 1st Qu.: -2.5773
Median : 5.69650 Median : 24.71 Median : 5.00 Median : -0.5739
Mean : 7.01900 Mean : 31.82 Mean : 7.34 Mean : 0.3210
3rd Qu.: 9.60323 3rd Qu.: 38.94 3rd Qu.:10.30 3rd Qu.: 2.6853
Max. :20.22214 Max. :107.38 Max. :36.00 Max. : 30.4535
zscore fold
Min. :-2.50740 Min. : 1.00
1st Qu.:-0.46053 1st Qu.: 8.25
Median :-0.11106 Median :15.50
Mean : 0.02631 Mean :15.50
3rd Qu.: 0.45814 3rd Qu.:22.75
Max. : 4.65631 Max. :30.00
plot(cv.exp$var1.pred ~ cv.exp$observed, cex = 1.2, lwd = 2)
abline(0, 1, col = "lightgrey", lwd = 2)
lm_exp <- lm(cv.exp$var1.pred ~ cv.exp$observed)
abline(lm_exp, col = "red", lwd = 2)r2_exp = summary(lm_exp)$r.squared
rmse_exp = hydroGOF::rmse(cv.exp$var1.pred, cv.exp$observed)
## 3) Modelo Gaussiano
cv.gau <- krige.cv(acumulado_24h ~ x + y, locations = pluvio, model = gau.fit)
summary(cv.gau)Object of class SpatialPointsDataFrame
Coordinates:
min max
x 634918.7 687966
y 7450222.8 7475508
Is projected: NA
proj4string : [NA]
Number of points: 30
Data attributes:
var1.pred var1.var observed residual
Min. : 0.471 Min. : 12.83 Min. : 0.00 Min. :-15.95270
1st Qu.: 2.869 1st Qu.: 14.63 1st Qu.: 2.20 1st Qu.: -2.33169
Median : 6.791 Median : 19.70 Median : 5.00 Median : -0.46101
Mean : 7.297 Mean : 29.96 Mean : 7.34 Mean : 0.04298
3rd Qu.: 9.971 3rd Qu.: 40.45 3rd Qu.:10.30 3rd Qu.: 1.88845
Max. :21.153 Max. :105.30 Max. :36.00 Max. : 30.93771
zscore fold
Min. :-2.765848 Min. : 1.00
1st Qu.:-0.464753 1st Qu.: 8.25
Median :-0.101036 Median :15.50
Mean : 0.004199 Mean :15.50
3rd Qu.: 0.248281 3rd Qu.:22.75
Max. : 4.660612 Max. :30.00
plot(cv.gau$var1.pred ~ cv.gau$observed, cex = 1.2, lwd = 2)
abline(0, 1, col = "lightgrey", lwd = 2)
lm_gau <- lm(cv.gau$var1.pred ~ cv.gau$observed)
abline(lm_gau, col = "red", lwd = 2)r2_gau = summary(lm_gau)$r.squared
rmse_gau = hydroGOF::rmse(cv.gau$var1.pred, cv.gau$observed)
## 4) Modelo Wave
cv.wav <- krige.cv(acumulado_24h ~ x + y, locations = pluvio, model = wav.fit)
summary(cv.wav)Object of class SpatialPointsDataFrame
Coordinates:
min max
x 634918.7 687966
y 7450222.8 7475508
Is projected: NA
proj4string : [NA]
Number of points: 30
Data attributes:
var1.pred var1.var observed residual
Min. : 1.979 Min. :14.15 Min. : 0.00 Min. :-11.4223
1st Qu.: 4.483 1st Qu.:20.05 1st Qu.: 2.20 1st Qu.: -4.7160
Median : 6.881 Median :27.94 Median : 5.00 Median : -0.5548
Mean : 7.024 Mean :28.33 Mean : 7.34 Mean : 0.3161
3rd Qu.: 9.391 3rd Qu.:35.76 3rd Qu.:10.30 3rd Qu.: 2.8856
Max. :13.217 Max. :48.97 Max. :36.00 Max. : 31.7935
zscore fold
Min. :-1.94021 Min. : 1.00
1st Qu.:-0.93925 1st Qu.: 8.25
Median :-0.13349 Median :15.50
Mean : 0.03151 Mean :15.50
3rd Qu.: 0.52837 3rd Qu.:22.75
Max. : 5.09689 Max. :30.00
plot(cv.wav$var1.pred ~ cv.wav$observed, cex = 1.2, lwd = 2)
abline(0, 1, col = "lightgrey", lwd = 2)
lm_wav <- lm(cv.wav$var1.pred ~ cv.wav$observed)
abline(lm_wav, col = "red", lwd = 2)# Criando uma tabela das estatística de R2 e RMSE
df.r2 <- data.frame(r2_lin, r2_exp, r2_gau, r2_wav)
df.rmse <- data.frame(rmse_lin, rmse_exp, rmse_gau, rmse_wav)
tabela <- data.frame(cbind(t(df.r2), t(df.rmse)))
colnames(tabela) <- c("R2", "RMSE")
rnames <- gsub("r2_", "", rownames(tabela)) # remove o prefixo r2 dos nomes das linhas
rownames(tabela) <- rnames # substitui o nome das linhas simplificadas na tab original
tabela R2 RMSE
lin 0.1482468633 6.882779
exp 0.0765441167 7.496573
gau 0.0680566292 7.723425
wav 0.0001772582 7.991183
Criando os grids do contorno da cidade do Rio de Janeiro para intermpolação
# Importando o contorno do Rio
contorno.rio <- shapefile("dados/chuva_rio/MUNIC_2K_2022_IPP_SIRGAS.shp")# Criando grade para interpolacao com a resolucao de 50m
r <- raster(contorno.rio, res = 50)
# Criando um objeto formato raster
rp <- rasterize(contorno.rio, r, 0)
# Trasnsformando o objeto raster no formato SpatialPixelsDataFrame
grid <- as(rp, "SpatialPixelsDataFrame")
sp::plot(grid)Krigagem
# Colocando os dados de chuva e o grid na mesma projecao
sp::proj4string(pluvio) = CRS(proj4string(contorno.rio))
mapa_chuva_lin <- krige(acumulado_24h ~ 1, pluvio, grid, model = lin.fit)[using ordinary kriging]
[using ordinary kriging]
[using ordinary kriging]
[using ordinary kriging]
Auto Krige
[using universal kriging]
krige_output:
Object of class SpatialPixelsDataFrame
Coordinates:
min max
x 623533.4 695333.4
y 7446574.5 7483374.5
Is projected: TRUE
proj4string :
[+proj=utm +zone=23 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_defs]
Number of points: 481750
Grid attributes:
cellcentre.offset cellsize cells.dim
s1 623558.4 50 1436
s2 7446599.5 50 736
Data attributes:
var1.pred var1.var var1.stdev
Min. :-0.01083 Min. : 0.07308 Min. : 0.2703
1st Qu.: 1.50803 1st Qu.: 14.16177 1st Qu.: 3.7632
Median : 5.25560 Median : 22.76226 Median : 4.7710
Mean : 6.66402 Mean : 24.90904 Mean : 4.7804
3rd Qu.:11.20533 3rd Qu.: 31.22818 3rd Qu.: 5.5882
Max. :35.96625 Max. :148.38326 Max. :12.1813
exp_var:
np dist gamma dir.hor dir.ver id
1 8 2254.389 11.19978 0 0 var1
2 27 4384.026 13.98185 0 0 var1
3 35 6073.421 36.23544 0 0 var1
4 54 8869.824 53.69680 0 0 var1
5 50 12016.434 41.83955 0 0 var1
6 47 15019.442 52.23264 0 0 var1
7 52 18659.164 85.79906 0 0 var1
var_model:
model psill range
1 Nug 0.0000 0.00
2 Exp 176.6964 33077.08
Sums of squares betw. var. model and sample var.[1] 0.0003263469
# Validação cruzada
auto.krige.cv <- autoKrige.cv(acumulado_24h ~ x + y, pluvio, model = "Exp")
summary(auto.krige.cv) [,1]
mean_error 0.36
me_mean 0.04904
MAE 4.544
MSE 57.79
MSNE 1.789
cor_obspred 0.2754
cor_predres -0.3772
RMSE 7.602
RMSE_sd 1.022
URMSE 7.594
iqr 4.636
Convertendo para o formato raster - auto krige
- Caso queira salvar a imagem raster em um arquivo formato geotiff para ler em algum SIG por exemplo:
Fazendo o mapa interativo com as estações
# Importando os dados referente as estatções pluviométricas
estacoes.sf <- read_sf("dados/chuva_rio/Estac_C3_B5es_Alerta_Rio.shp")
# Convertendo UTM para Lat Long das estacoes
estacoes.longlat <- st_transform(estacoes.sf, "+proj=longlat +ellps=WGS84 +datum=WGS84")
estacoes.longlat$coords <- st_coordinates(estacoes.longlat)
estacoes.longlat$X <- estacoes.longlat$coords[, 1]
estacoes.longlat$Y <- estacoes.longlat$coords[, 2]# Importando a malha de bairros
bairros.sf <- read_sf("dados/chuva_rio/BAIRROS_2K_2022_IPP_SIRGAS.shp")
# Convertendo UTM para Lat Long a malha dos bairros
bairros.longlat <- st_transform(bairros.sf, "+proj=longlat +ellps=WGS84 +datum=WGS84")# Convertendo o raster da chuva para lat long
raster_chuva_longlat <- projectRaster(raster_chuva, crs = CRS("+proj=longlat +datum=WGS84"))# Definindo Paleta de cores da superfície interpolada da chuva
pal <- colorNumeric(c("#000066", "#00c8f8", "#F0E68C", "#FFFF00", "#FF8C00"), values(raster_chuva_longlat),
na.color = "transparent", reverse = T)# Construindo o mapa interativo via leaflet
leaflet(data = estacoes.longlat, options = leafletOptions(attributionControl = FALSE)) |>
# addTiles() |>
addProviderTiles("CartoDB.Positron", group = "Ruas") |>
addProviderTiles("Esri.WorldImagery", options = providerTileOptions(opacity = 0.7),
group = "Satélite") |>
addProviderTiles(providers$CartoDB.Voyager, group = "Voyager") |>
addProviderTiles(providers$Stamen.Toner, group = "Toner") |>
setView(lng = -43.42, lat = -22.9, zoom = 10.4) |>
# addProviderTiles(providers$CartoDB.Voyager) |>
addMarkers(~X, ~Y, popup = ~as.character(est), label = ~as.character(est), group = "Estações") |>
############## Polígonos dos Bairros ################
addPolygons(data = bairros.longlat, weight = 3, color = "darkblue", smoothFactor = 1,
fill = FALSE, labelOptions = labelOptions(style = list(`font-weight` = "normal",
padding = "3px 8px"), textsize = "13px", direction = "auto"), group = "Bairros") |>
########## Adicionando o raster #########################
addRasterImage(raster_chuva_longlat, colors = pal, opacity = 0.8, group = "Chuva: 1 semana") %>%
leaflet::addLegend(pal = pal, values = values(raster_chuva_longlat), title = "Chuva Acumulada - 1 semana",
group = "Chuva: 1 semana") |>
############## Controle das layers (botoes) ################
addLayersControl(baseGroups = c("Voyager", "Ruas", "Satélite", "Toner"), overlayGroups = c("Estações",
"Bairros", "Chuva: 1 semana"), options = layersControlOptions(collapsed = FALSE),
position = "bottomleft") |>
########## Desabilitando os grupos ################
hideGroup(group = c("Bairros"))Bibliografia sugerida
Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds). Análise Espacial de Dados Geográficos. Brasília, EMBRAPA, 2004.
Interactive Spatial Data Analysis by Trevor C. Bailey , Anthony C. Gatrell Routledge, 1995
Applied Spatial Statistics for Public Health Data; Lance A. Waller, Carol A. Gotway Wiley-Interscience 1St ed. 2004
Applied Spatial Data Analysis with R; Roger S. Bivand, Edzer Pebesma , Virgilio Gomez-Rubio Springer; Edição: 2nd ed. 2013
Geocomputation with R by Robin Lovelace, Jakub Nowosad and Jannes Muenchow. Geocomputation with R, 2021.
Spatial Statistics Workbook of Department of Criminology at the University of Manchester by Reka Solymosi and Juanjo Medina Crime Mapping in R
Spatial Data Science with R Spatial Data Science with R
GeoComputation and Spatial Analysis practicals by Lex Comber On-line Book
GEOG5917 Big Data & Consumer Analytics - RStudio Practicals by Lex Comber On-line Book